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We study the dynamics of phase synchronization in growing populations of discrete phase oscil¬ 
latory systems when the division process is coupled to the distribution of oscillator phases. Using 
mean field theory, linear stability analysis, and numerical simulations, we demonstrate that cou¬ 
pling between population growth and synchrony can lead to a wide range of dynamical behavior, 
including extinction of synchronized oscillations, the emergence of asynchronous states with unequal 
state (phase) distributions, bistability between oscillatory and asynchronous states or between two 
asynchronous states, a switch between continuous (supercritical) and discontinuous (subcritical) 
transitions, and modulation of the frequency of bulk oscillations. 


I. INTRODUCTION 

Synchronization phenomena in collections of coupled 
oscillators and excitable elements are widely studied in 
statistical physics [H-Q, in part because they represent 
prototypical nonequilibrium phase transitions exhibiting 
time-translational symmetry breaking. In addition to 
their theoretical value, models of synchronization offer in¬ 
sight into a diverse collection of physical, chemical, and 
biological phenomena ranging from bulk oscilla¬ 

tions in chemical reactions to phenotypic or behavioral 
synchronization in populations of living organisms. 

Synchronization plays a particularly important role in 
biological systems, where coherent oscillations may serve 
biological or behavioral function. Exaimles abound, in¬ 
cluding rhythmic flashing of fireflies [3|, coherent be¬ 
havior of neurons in human neocortex Q or primate 
retina Q, circadian oscillations in cyanobacteria 
vertebrates 11211 and mammals [l3lll4|. an d sy nchronized 
cell division [l5j or protein dynamics [1^ llZl in popula¬ 
tions of single-cell organisms. In many of these systems, 
the timescales of the oscillations are well-separated from 
the timescales of population growth, allowing one to ne¬ 
glect its effects on macroscopic synchrony. In turn, the 
vast majority of theoretical studies on coupled oscilla¬ 
tors have dealt with a fixed population size. However, 
this restriction may not always be applicable, as an in¬ 
creasing number of systems have been shown to exhibit 
oscillations that occur on similar timescales as popula¬ 
tion growth [Ml EMI. This overlap in timescales 
raises the question of how population growth might af¬ 
fect synchronization, given that the population wide dis¬ 
tribution of oscillator phases may be strongly coupled to 
the growth process. 

In this work, we explore the effects of population 
growth and the corresponding redistribution of phases 
on the synchronization properties of a simple class of 
models for both coupled oscillators and excitable ele¬ 
ments. While the paradigmatic Kuramoto model has 
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paved the way to much of our current understanding of 
synchrony ii, oscillator models with discrete phases 
have gained increasing attention because of their relative 
mathematical simplicity [20l - [^ . Here, we focus on these 
discrete phase models because they can be readily modi¬ 
fied to integrate oscillator growth and the potential phase 
dependence of the division and birth processes. 

Using a combination of numerical simulations, mean 
field theory, and linear stability analysis, we find that the 
redistribution of phases induced by population growth 
can disrupt synchronization via either continuous (super¬ 
critical) or discontinuous (subcritical) transitions in both 
discrete phase oscillators and excitable elements. We ob¬ 
serve a range of dynamical behaviors, including bistabil¬ 
ity between two asynchronous states or between asyn¬ 
chronous and oscillatory states, a switch between super¬ 
critical and subcritical transitions as growth is increased, 
the existence of asynchronous states with unequal phase 
distributions, or modulation of the bulk oscillation fre¬ 
quency. These results demonstrate that even in minimal 
models, the coupling between population growth and os¬ 
cillator phase can profoundly affect synchronization and 
even lead to new dynamical states that do not exist in 
the absence of this coupling. 

The paper is divided into two sections, with the first 
devoted to discrete phase oscillators and the second to 
excitable, discrete phase systems. In Section III A1 we 
briefly review the discrete phase oscillator model and ex¬ 
tend it to capture population growth. In Section lllBI we 
develop a mean field approximation, and in sections III Cl 
and ED] we use linear stability analysis and numerical 
simulations to explore the effects of growth when divi¬ 
sion is independent of state or strongly state-dependent, 
respectively. In Section llH Al we outline a model for grow¬ 
ing populations of excitable elements, in Section llH BI we 
describe numerical simulations of the model, and in Sec¬ 
tion IIII Cl we use mean field theory and linear stability 
analysis to derive complete phase diagrams for the grow¬ 
ing populations. We conclude with a discussion of the 
results in Section HVl 
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II. DISCRETE PHASE OSCILLATORS 
A. Model for Growing Oscillator Populations 

We will model an active oscillator as an m-state system 
governed by unidirectional transitions between states, 
1 —>• 2 —>• 3... —>• m —^ 1 —^ 2... [ 2 ^. The states repre¬ 

sent a type of discretized phase, and formally, the state 
space of a single oscillator can be described by a phase 
variable 0 < (^ < 27r in the limit m —> 00 . For sim¬ 
plicity, we restrict ourselves here to finite m and study 
a discrete phase oscillator with the minimum number of 
states (m = 3) required to generate a Hopf bifurcation 
and, hence, macroscopic oscillations in a coupled popu¬ 
lation [2l|. In addition to their relative simplicity, these 
discrete phase models are often appropriate descriptions 
for biological or chemical systems, where the oscillations 
commonly occur on a discrete state space. 

In the absence of coupling, each oscillator transitions 
irreversibly between states (1 —)> 2 —>■ 3 —)> 1) with a 
probability per unit time g, which sets the oscillator’s in¬ 
trinsic frequency. The model is therefore an example of a 
Markov chain. Coupling between oscillators is achieved 
by allowing these transition rates to depend on the frac¬ 
tion of the oscillators in each state; hence, g is replaced 
by a function T^, which is the probability per unit time 
for a given oscillator to transition from state * — 1 to 
state i (with i = 1,2,3 (mod 3)). T^ also depends on 
a real parameter a > 0, which measures the strength 
of the coupling between oscillators (hence dVi/da > 0). 
While a number of nonlinear coupling functions have 
been used in previous studies pH . Wi \ , for now we will 
leave the coupling function unspecified but explicitly note 
its dependence on the fraction of oscillators Pi in state i, 
hi = Ti{Pi). The primary requirement for this function 
is that it facilitates phase coherence between oscillators 
and leads to a Hopf bifurcation at a positive value of 
a. We discuss these constraints in more detail in what 
follows. 

In the mean field limit of all-to-all coupling, the frac¬ 
tion of oscillators Pi in state i is governed by the contin¬ 
uous time master equation. 

Pi = ~Pi^i+l + Pi-l^i- (1) 

The total number of oscillators is conserved in this model 
= 1), and the fixed point P* = 1/3 becomes 
unstable via a Hopf bifurcation as long as 


above some critical value of a = Oc P3 • Equation [2l 
r = r(x,a )|^^^/3 and T' = For larger a > 

ttc, macroscopic oscillations occur. In what follows, we 
consider the class of oscillator models where Equation [2] 
is satisfied for some a > Uc > 0. 

To incorporate population growth, we introduce two 
minimal mechanisms by which growth and oscillator 


phase may be coupled. First, we allow each oscillator in 
state i to give birth to a new oscillator with probability 
per unit time tik, with ti = m = Z. The dimension¬ 
less weighting factors couple the rate of division to the 
state of the oscillator, and fc is a growth rate constant. 
In this work, we restrict our attention to two limiting 
cases: phase independent growth, = 1, and strongly 
phase-dependent growth, = 3(5^, where 5ij is the Kro- 
necker delta. In the former case (Case I), oscillators in 
all states are equally likely to divide, so division itself is 
independent of the state of the oscillator. In the latter 
case (Case 2), division can only occur for oscillators in 
state j. We choose j = 1 without loss of generality. This 
state dependence of division could be relevant in a num¬ 
ber of biological settings, such as synchronization in cell 
division itself, where the oscillators in certain states-for 
example, those in certain stages of the cell cycle-divide 
preferentially. State-dependence of division could also 
occur in oscillations of protein dynamics, as transcrip¬ 
tional and translational activity are often linked to the 
cell cycle. 

Second, we allow the phase of the daughter oscillator 
to depend probabilistically on the state of the mother. 
Specifically, with each division, the new daughter oscil¬ 
lator has probability y to be in the same state as the 
mother and a probability (1 — x)/2 to be in each of the 
two remaining states. This assumption could again be 
relevant in many biological applications, where cell divi¬ 
sion can lead to a repartitioning of intracellular contents, 
such as proteins, that may be fundamental to the oscil¬ 
lation. We stress that our goal here is not to incorpo¬ 
rate biological details into a system-specific model, but 
rather to introduce state-dependent division and birth in 
a minimal model. While other rules for coupling single 
oscillator dynamics to division are possible, we will show 
below that the simple rules above lead to a rich collection 
of dynamical behaviors. 


B. Mean Field Theory 

To develop a mean field approximation for this model, 
we begin by writing evolution equations for W, the num¬ 
ber of oscillators in state i, as 

Ni = —WEi+i + Ni_iri+ 

k (3) 

ukNiX + — (ei+iAi_|_i -I- ei_iW-i)(l ~ x) 

with i = 1,2,3 (mod 3). The first two terms capture 
the nonlinear coupling between oscillators that drives 
synchronization, analogous to the two terms in Equa¬ 
tion [1] while the latter two terms account oscillator divi¬ 
sion, as described in Section Hi Al Since the total number 
of oscillators = -^(0) rewrite Equation |3] in 

terms of Pi, the fraction of oscillators in state i, using 
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Pi = Ni/N — NiN/N'^, to arrive at 

Pi = —Pi^i+l + Pi-l^i + iikPiXP 

k N (4) 

^2 (^i+l-fj+l P £i—l-Pi-l)(l ~ x) ~ 

where N/N = follows directly from Equa¬ 

tion [S] The model is now fully specified by two differen¬ 
tial equations (e.g. for Pi and P^) and the normalization 
condition P 1 -I-P 2 + T 3 = 1- Without loss of generality, we 
can set 5 = 1 , which is equivalent to measuring time in 
units of g and replacing k with k/g] in what follows, we 
use k for economy of notation. Equation|T]provides a gen¬ 
eral mean field description valid in the limit of all-to-all 
coupling and can be solved numerically for any choice of 
parameters. To make further analytical progress, we now 
restrict our attention to the limiting cases of phase inde¬ 
pendent growth, Ci = 1, and strongly phase-dependent 
growth, Cj = S6ij. 


C. Case 1: Division is independent of state 


When division takes place independently of the state 
of each oscillator (e^ = 1), Equation 0] reduces to 


Pi — —Pi^i+l + Pi-l^i — kPi(l — x)P 

k ( 5 ) 

Several things become apparent from inspection. Eirst, 
for X = 1 (daughter oscillators are always in the same 
state as the mother). Equation [5] reduces to Equation [T] 
In this case, oscillator division leads to an exponentially 
increasing number of oscillators over time, but all syn¬ 
chronization properties-which depend on the fraction of 
oscillators in each state, not the total number-will remain 
unchanged. Second, it is clear that the asynchronous 
fixed point Pi = 1/3 is a solution to Equation [5] for all 
parameter values. 

To analytically explore the effects of division on os¬ 
cillator synchrony, we linearize around the asynchronous 
fixed point; the corresponding Jacobian matrix, J, eval¬ 
uated at the fixed point, is 


J = 


-|(i-x)fc-2r + ir' 

r + ir 


-T - f-T' 


-§(i-x)A:-r + |r 


( 6 ) 

where T and T' are the coupling function and its deriva¬ 
tive, respectively, evaluated at Pi = 1/3. The matrix J 
has complex conjugate eigenvalues X = aPiuj, with 


« = ^(-3(i-x)fc-3r + r') 


CU 




=r' 


which cross the imaginary axis at 

1 


k = kc = 


3(1 -X) 


V3 

is at 

(-3r(a,) + r'(a,)), 


( 7 ) 


( 8 ) 
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FIG. 1. Top panels: values of successive maxima / minima 
of the time series Pi(t) (black, dark) and P 2 (t) (red, light) 
for each value of the dimensionless growth k for model in 
Equation [S] Top left, x = Oi bottom left, x = 1/3; top 
right, X = 2/3; bottom right, x = 1- Bottom panel: order 
parameter q for x = 0 (light blue), x = 1/3 (black), x = 2/3 
(dark blue), and X = 1 (red). Different shapes represent 
different intrinsic oscillator frequencies {g = 0.75, circles; g = 
1, crosses; g = 1.5, triangles). In all panels, a = 20 > a^.o. 
No — 500, initial conditions are given by Pi = 2/3, P 2 = 1/3, 
and the coupling function F; = 1 -|- oPf. 



where we have explicitly noted the dependence of T 
and r' on the critical value of the coupling, Oc. Equa¬ 
tion [5] provides a relationship between the critical val¬ 
ues of growth kc and coupling Oc, which separate syn¬ 
chronous and asynchronous behavior. For A: = 0, Equa¬ 
tion [5] reduces to Equation 01 More generally, for k > kc, 
the fixed point is stable and no oscillations occur; suffi¬ 
ciently fast population growth therefore disrupts other¬ 
wise synchronized populations. In addition, at the bifur¬ 
cation point, the frequency of the oscillations is given by 
ui = luq = ^ ^\/3r(ac) + ■^r^(®c)^ and is therefore ex¬ 
pected to be modified by the addition of growth, which 
changes Qc- 
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FIG. 2. Top panels: values of successive maxima / minima 
of the time series Pi{t) (black, dark) and P 2 {t) (red, light) 
for each value of the dimensionless growth rate k for model 
in Equation [H Top left, X = 0; bottom left, x = 1/3; top 
right, X = 2/3; bottom right, X = 1- Bottom panel: order 
parameter q for x = 0 (light blue), x = 1/3 (black), x = 2/3 
(dark blue), and X = 1 (red). In all panels, a = 4 > Oc.o, 
No = 500, g = 1, initial conditions are given by Pi = 2/3, 
P 2 = 1/3, and the coupling function Vi = 


To determine the nature of the bifurcation (subcritical 
or supercritical), we calculate the first Lyapunov coeffi¬ 
cient, li- The sign of /i, which is analogous to the coef¬ 
ficient of the third order term in the normal form for a 
Hopf bifurcation, is negative for supercritical and positive 
for subcritical Hopf bifurcations. For an n-dimensional 
dynamical system x = f{x, e) with an equilibrium point 
X = x^ undergoing a Hopf bifurcation at parameter value 
e = , li can be calculated following the projection pro¬ 

cedure in [33 as 

h =^Re[{p, C{q, q, q)) - 2(p, B{q, A-^B{q, q))) 

ZujQ ( 9 ) 

+ (p, B{q, {2iuJoI - A)-^B{q, g)))], 
where (.,.) is the complex scalar product, / is the identity 
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FIG. 3. Upper panel: frequency of macroscopic oscillations 
for model in Equation [5l corresponding to the maximum peak 
of the power spectrum of Pi{t) for x = 0,1/3, 2/3,1 (upper 
triangles, circles, and squares, respectively). Curves originat¬ 
ing at tj « 6.5 correspond to Fi = curves originating 

at Lo ~ 5.2 correspond to Fi = 1 -|- aP/. Lower panel: order 
parameter, q, as a function of kc — k near the critical point. 
Open shapes correspond to Fi = 1 -|- aPf ; small closed circles 
correspond to Fi = Open shapes: X = 0 (light blue), 

X = 1/3 (black), x = 2/3 (dark blue) at different intrinsic os¬ 
cillator frequencies [g = 0.75, circles; <7 = 1, crosses; g = 1.5, 
triangles). Closed circles: g = 1 and x = 0 (blue), x = 1/3 
(red), X = 2/3 (green). In all panels. No = 500 and initial 
conditions are given by Pi = 2/3, P 2 = 1/3. o = 20 > flc.o 
for Fi = 1 -|- aP/; a = 4 > Uc.o for Fi = Dashed line: 

mean field scaling q ~ {kc — k^A^ Curves are shifted slightly 
to allow for visualization of all curves simultaneously. 


matrix, and p and q are the right and left eigenvectors, 
respectively, of the Jacobian A = given by 


Aq = iujoq, 

A^p = —icoop. 


( 10 ) 


The vectors are normalized so that {p,q) = 1, and 
the functions B(u,v) and C{u,v,w) are multilinear, n- 
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dimensional vector functions given by 


B{u,v) = ^ 


3,k=l 


dtpjdipk 


UjVk 




C{u,v,w)= ^ 






( 11 ) 


UjVkWl. 




Specifically, for the model in Equation!^ we find 
'-l-iy/S 1 \ 


9 = 


p = \ -i 


2 V 2 ' V2J’ 

(2 3V2 + zV6^ 

3’ 6 


( 12 ) 


independent of the coupling function E. Following 
straightforward algebraic manipulations, we arrive at 


h = 


_ V3 fr"' - er'' 


4 V 3r + r' 


(13) 


where primes indicate derivatives of the coupling function 

r = r(a:,ac)|^^,/ 3 . 


^ dr{x,ac) I 

dx \x=i/v 

d‘^V{x,ac) I 

Q^2 |£c= 1 / 3 ’ 

i9^r(a:,ae) , 

dx^ la;=l /3 


r" = 

p/// _ 


(14) 


Interestingly, Equation [T3] suggests that the nature of the 
Hopf bifurcation is determined by the magnitude of the 
derivatives of the coupling function at the critical point. 
Because increasing growth rate, k, will change the critical 
value Oc, it is possible for growth to not only change the 
coupling required for synchronization, but also the na¬ 
ture of the transition itself. For example, if we make the 
physically realistic assumptions that E^ > 0 and E' > 0, 
the condition E'" = 6E" separates supercritical and sub- 
critical transitions. To illustrate this point, in the next 
section we consider two specific examples of the coupling 
function and show that both supercritical and subcritical 
bifurcations are possible, depending on its derivatives. 


1. Examples of supercritical and subcritical growth-induced 
bifurcations 

In this section, we study two specific coupling functions 
to demonstrate the rich dynamics possible in this class 
of models. First, we consider a coupling of the form 

riP,,a) = l + aP^, (15) 

which satisfies Equation [5] for a > Oc = 9 in the absence 
of population growth. For nonzero growth rate fc, we 



FIG. 4. Order parameter g as a function of k for Equation [5] 
with Ei = e“^* and a = 6.75. Simulations were run from 
two sets of initial conditions: Pi = 2/3, P 2 = 1/3 (red) and 
Pi = 1/3, P 2 = 1/3 (blue), leading to different steady state 
behavior. Insets: phase portraits for A: = 12 for each set of 
initial conditions, x = 0 and No = 5000 for all points. Dashed 
lines indicate region of bistability. 


can rearrange Equation [5] to show the critical value Oc is 
increased to 


flc = Oc.o + 9fc(l - x), (16) 


where ac,o is the critical coupling in the absence of 
growth. In a synchronized population (a > Cc^o), intro¬ 
ducing population growth therefore quashes the macro¬ 
scopic oscillations when 


k > kc 


a — 9 

9(1-x)' 


(17) 


At the transition point, the frequency of oscillations is 
given by 



(ttc + 3) = \/3 



^fcc(l-x)^ ■ 


(18) 


Hence, for a fixed value of fc = kc, the frequency at 
the transition point decreases linearly with x, eventu¬ 
ally reaching the transition frequency ujq = 2\/3 of the 
non-dividing model. On the other hand, for a fixed value 
of a > ttcfi, the frequency at the transition point ap¬ 
proaches the value uj = ^(a -|- 3), independent of X) as 
k approaches kc- 

Finally, we can calculate the first Lyapunov coefficient, 
h, which depends on fc, as 


9^(l + fc(l-x)) 
4-f 3A:(l-x) 


(19) 
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Hence, /i < 0 for all k and 0 < x < 1, indicating that 
the bifurcation remains supercritical in the presence of 
population growth. 

As a second example, we consider the coupling 

r(P„a) = e“^-. (20) 

In the absence of growth, the model undergoes a Hopf 
bifurcation at a = ac,o = 3. In the presence of growth, 
the critical value Oc increases and is given by the solution 
to 


3fc(l-x) = (a-3)e“/^ (21) 


Equivalently, if we start with a synchronized population 
at a given value of a > Oc^O) the oscillations will be de¬ 
stroyed when k > kc = {a — 3)e“/^/(3(l — x))- The 
frequency at the transition point is given by 


(uc -b 3)e“/3 

27 ! 


( 22 ) 


which increases when growth is introduced and Uc > ac,o- 
Most interestingly, the first Lyapunov coefficient, /i, is 


\/3(Qe - &)al 
4(3 -\- Qc) 


(23) 


which changes sign at Oc = 6 (or equivalently, when V" = 
6r'). For small growth rates such that ac,o < Oc < 6, the 
bifurcation between synchronous and asynchronous state 
is supercritical. However, for larger growth rates, Uc > 6 
and the bifurcation becomes a subcritical, discontinu¬ 
ous transition. In this case, population growth-and the 
corresponding redistribution of oscillator phases-leads to 
a fundamentally different transition for sufficiently high 
growth rates. 

To confirm these results numerically, we simulated 
growing populations of oscillators using the Gillespie al¬ 
gorithm [33| for a wide range of k (growth rate), g (oscil¬ 
lator natural frequency), and x (probability of daughter 
being in same state as mother). Unless otherwise noted, 
simulations were run starting with Nq = 500 oscillators 
that are strongly coupled (a > ac,o)- For nonzero fc, 
the total number of oscillators grows approximately ex¬ 
ponentially, so computer memory limits simulations to 
relatively short time periods. To circumvent this lim¬ 
itation, we allowed simulations to run until the total 
number of oscillators reached Nmax = lOAb; when the 
number of oscillators exceeded Nmax, we automatically 
reset the total number of oscillators to Nq while preserv¬ 
ing the fractional distribution of oscillator states. While 
this numerical procedure effectively underestimates the 
fluctuations observed in the oscillations, our goal is to 
approximate a thermodynamic limit N ^ oo^ and even 
the modest starting number N = 500 yields relatively 
small fluctuations in macroscopic behavior. 

After simulations have reached steady state, we visu¬ 
alize the dynamics by plotting the values of successive 


maxima / minima of the time series Pi{t) (black, dark) 
and P 2 it) (red, light) for each value of the dimensionless 
growth rate k (Figures [T] and [2]). We also calculated the 

synchronization order parameter q = ((|Z — {Z)t\'^)t) , 

which was originally proposed in ( 2 ^ . In this definition, 
— W i® tFe (not yet averaged) Kuramoto 

order parameter [1 i, Oj = 27rk/3, k = ( 0 , 1 , 2 ) is the 
discretized phase of oscillator j, and angle brackets repre¬ 
sent an average over time. The standard Kuramoto order 
parameter is not appropriate for models where rotational 
symmetry is absent and, consequently, non-oscillating 
steady states can lead to nonzero values of the order pa¬ 
rameter, despite the absence of collective synchrony. As 
shown in [ 22 j |. the order parameter q is akin to a gener¬ 
alized standard deviation of Z{t) and removes the bias 
due to a lack of rotational symmetry. Nonzero values of 
q correspond to synchronized, oscillatory states. 

For Ti = 1 -b aPf at a = 20 > ac,o (Figured]) and 
Vi = at o = 4 > ac,o (Figure d]), oscillations 
smoothly decrease in amplitude with increasing k, lead¬ 
ing eventually to a completely asynchronous state. The 
value of the critical growth k is consistent with the linear 
stability analysis, Equation dZI and Equation dT] it in¬ 
creases with X until, at x = Ij the oscillations are undis¬ 
turbed by even large growth rates. 

Our numerical simulations indicate that population 
growth affects not only global synchronization, but also 
the frequency of the macroscopic oscillations in the syn¬ 
chronized state. To explore this frequency dependence 
systematically, we calculated the power spectrum of the 
time series Pi{t) for each simulation in steady state. Fig¬ 
ure [3] (top panel) shows the frequency of the dominant 
peak in the power spectrum, which characterizes the fre¬ 
quency of the macroscopic oscillations. In both examples, 
increasing k leads to a monotonically increasing oscilla¬ 
tion frequency until, at sufficiently high values of k, the 
macroscopic oscillations are no longer present. 

Because these models are globally coupled and exhibit 
supercritical bifurcations, one expects that the order pa¬ 
rameter should scale as q ^ {kc — k)^ near the critical 
point, with /3 = 1/2 the standard mean field scaling ex¬ 
ponent. Indeed, we find that in our numerical simula¬ 
tions, the order parameter approximately follows mean 
field scaling near the critical point, independent of the 
value of X or the specific coupling function chosen (Fig¬ 
ure d] bottom panel; note that we have slightly shifted 
the curves relative to one another so that each is visi¬ 
ble). One would expect similar behavior for any choice 
of Ti where a supercritical transition occurs as k eclipses 
a critical value kc- 

By contrast, the transition becomes discontinuous 
when Ti = and a is increased beyond a = 6. 

For example, at a = 6.75, growth induces a subcrit¬ 
ical Hopf bifurcation (Figure dj), as indicated by the 
discontinuous drop in the order parameter and the cor¬ 
responding bistability. Interestingly, in the bistable re¬ 
gion, which we find to exist for 11.9 < k < 12.1, syn¬ 
chronous oscillations co-exist with an asynchronous fixed 
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point. For populations initiated with approximately uni¬ 
form phase distributions, the system settles into a stable 
asynchronous state; for highly non-uniform distributions, 
the population undergoes stable oscillations indicative of 
synchrony (upper right inset). Similar behavior would be 
expected in the class of models where li can switch signs 
as growth increases (see Eauation lldl) . indicating that in 
this class of models, coupling between population growth 
and phase can modify the nature of the Hopf bifurcation. 




FIG. 5. Top panels: values of successive maxima / minima of 
the time series Pi(t) (black, dark) and 7*2(t) (red, light) for 
each value of the dimensionless ratio k for the model Equa¬ 
tion [211 Top left, X = 0; bottom left, x = 1/3; top right, 
X = 2/3; bottom right, x = 1- Bottom panel: order parame¬ 
ter q for X = 0 (light blue), x = 1/3 (black), x = 2/3 (dark 
blue), and X = 1 (red). Different shapes represent different in¬ 
trinsic oscillator frequencies {g = 0.75, circles; g = 1, crosses; 
g = 1.5, triangles). Inset: frequency ui of macroscopic oscilla¬ 
tions, corresponding to the maximum peak of the power spec¬ 
trum of Pi{t) for X = 0i 1/3, 2/3, 1 (upper triangles, circles, 
squares and stars, respectively). In all panels, a = 20 > a^.o, 
Nq = 500, initial conditions are given by Pi = 2/3, P 2 ~ 1/3, 
and the coupling function Fi = 1 -|- aP/. 


D. Case 2: Division occurs only in one state 

When division occurs only in one state (e^ = 3<5ii), 
equation [4] reduces to 

Pi = —PiPi+l + Pi-lPi + “ikPiidi^xX — Pl) + 

3 A- (24) 

y(<5,,3Pi+l+<5i,2Pz-l)(l-x) 

where Sij is the Kronecker delta. Because the choice of 
has broken the rotational symmetry of the model, the so¬ 
lution Pi = 1/3 is not, in general, a steady state solution 
to Equation [M] In what follows, we first consider the 
case X = 1/3, which is amenable to analytical treatment, 
and then go on to numerically explore the general case 
0<X< 1- 

When X = 1/3, Equation (Ml reduces to 

■ ^ (25) 


Pi — —Pi^i-\-l + Pi-l^i 


3kP,, ( 


for which Pi = 1/3 is always a solution. The correspond¬ 
ing matrix J is identical to Equation [5] when x = 1/3, 
and the linear stability analysis yields identical results to 
those in Equations IMT^ For example, if F^ = 1-|-aP^^, we 
have a critical coupling of Oc = Oc.o + 6 fc, and when the 
coupling is fixed at a = 20 , synchronization is destroyed 
when k > kc = 11 / 6 , which is consistent with numeri¬ 
cal simulations (Figure [5]). In addition, at the transition 
point w = ^(oc -b 3) = •\/3(2 -b kc) and the transition 
is always supercritical ^ ■ Hence, 

while Equation [5] and Equation [24] correspond to micro¬ 
scopically distinct mechanisms and differ in higher order 
terms, when x = 1/3 the linear stability properties of 
the fixed point Pi = 1/3, which determine synchroniza¬ 
tion properties near the phase transition, are identical. 
Intuitively, this correspondence arises because, in both 
cases, the division process redistributes the oscillators 
uniformly between the 3 possible states. 

For other values of x, however, state-dependent di¬ 
vision can give rise to entirely new dynamics. To ex¬ 
plore this behavior, we performed numerical simulations 
for a wide range of parameters, as in Section III Cl As 
a prototype model, we choose F^ = 1 -b aPi, but we 
later show that the similar behavior is observed for other 
coupling functions and even in continuous phase models. 
While the amplitude of the oscillations decreases with in¬ 
creasing growth rate k (Figure |S]), the oscillations of Pi 
and P 2 do not always occur around the symmetric val¬ 
ues (PijPz) = (1/3,1/3). Furthermore, for sufficiently 
large values of fc, the system appears to settle into a non¬ 
oscillating fixed point where the fraction of oscillators 
in states 1 and 2 can be significantly different. As an 
example, for x = 2/3 (Figure [SI upper right) the oscilla¬ 
tions cease at fc ~ 3, leading initially to a non-oscillating 
state where Pi < P 2 . This is somewhat counterintuitive, 




















FIG. 6. Upper left panel: Steady state solution P* of the 
mean field model, Equation 1241 Stability is indicated by 
marker (open circles, unstable points; closed circles, stable 
points; squares, saddle points). Dashed box is region where 
synchronous oscillations stably coexist with a non-oscillating 
state. Solid box is region where two asynchronous states sta¬ 
bly coexist. Upper right panel: Two example time series 
of P\ from numerical simulations with Nq = 2000 starting 
from identical initial conditions. Pi = 0.45, P 2 = 0.45. In 
one case, the system settles into stable oscillations (red). In 
the other case, a fluctuation drives the system to the asyn¬ 
chronous fixed point following initial oscillations (black). Bot¬ 
tom panel: phase portraits for k = 2.3, 2.6, 2.9, 3.2, 3.5 (left to 
right). Thin (blue) lines are example trajectories, thick (red) 
lines are stable limit cycles; stability of fixed points is indi¬ 
cated by markers (open (red) circles, unstable; closed (red) 
circles, stable; squares (black), saddle). In all panels, y = 1. 


as only oscillators in state 1 divide, and the majority of 
daughter cells (y = 2/3) are also in state 1. As fc is fur¬ 
ther increased, the population is eventually dominated by 
oscillators in state 1, as one might expect. Each individ¬ 
ual oscillator continues to cycle through all three states, 
but on average, the distribution of states is not uniform. 

In addition, growth can dramatically affect the fre¬ 
quency of oscillations in the synchronous state, but the 
effect is no longer monotonic for all values of y. As fc is 
increased, we find oscillations of increasing frequency for 
y = 0 and y = 1/3, approximately constant frequency 
for y = 2/3, and rapidly decreasing frequency for y = 1 
(Figure m lower panel inset). It is somewhat surprising 
that the choice of y, alone, can lead to either an increase 
or decrease in the overall oscillation frequency. Interest¬ 
ingly, the case y = I also appears to undergo an abrupt 
transition for k ~ 2.9, indicating that the transition is 
qualitatively different from the supercritical Hopf bifur¬ 
cation in the non-growing model. 

To examine this transition in detail, we performed lin- 




P2 


FIG. 7. Upper panel: Steady state solution P* of the 
mean field model. Equation 1241 Stability is indicated by 
marker (open circles, unstable points; closed circles, sta¬ 
ble points; squares, saddle points). Dashed box is region 
where synchronous oscillations stably coexist with a non¬ 
oscillating state. Solid box is region where two asynchronous 
states stably coexist. Bottom panel: phase portraits for 
k = 2.7, 2.85, 3.0, 3.15, 3.3 (left to right). Thin (blue) lines are 
example trajectories, thick (red) lines are stable limit cycles; 
stability of fixed points is indicated by markers (open (red) 
circles, unstable; closed (red) circles, stable; squares (black), 
saddle). In all panels, y = 1. 


O' 



FIG. 8. Order parameter q in growing population of identical 
Kuramoto oscillators, each with intrinsic frequency oj = 10. 
Goupling strength K = 1\ Number of initial oscillators Nq = 
500. When population size reaches N = 10®, the population 
is reset to a size of A = 1000 while preserving the distri¬ 
bution of oscillator phases. Points represent averages from 
simulations of 10 independent trials. Insets: Time series of 
order parameter Z{t) for two typical simulations {k — 0.32 
in both cases) in the bistable region; the population stochas¬ 
tically switches between two stable states, one synchronized 
and one asynchronous. 
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ear stability analysis numerically for y = 1. In particular, 
for y = 1 , we find that the fixed point P* remains close 
to Pi = 1 /3 for A: < 3 and loses stability via a supercrit¬ 
ical Hopf bifurcation at A: Ri 3, in apparent contradiction 
with numerical results, which indicate a discontinuous 
transition at a smaller value of k (see Figure [5|). Inter¬ 
estingly, a more thorough numerical analysis reveals that 
a second branch of stable solutions arises at A; r; 2.5 (Fig¬ 
ure upper left). The emergence of this solution branch 
leads to two types of novel bistable behavior and un¬ 
derlies the abrupt transition observed in numerical sim¬ 
ulations. As k is increased above k r; 2.5, one finds 
bistability between a synchronized oscillating state and 
a non-oscillating state with Pi > P 2 > P 3 . Surprisingly, 
the synchronized oscillations occur only when oscillator 
states are initially distributed within a relatively small 
section of phase space. Initial conditions with Pi = I for 
i =1, 2, or 3, for example, will (counterintuitively) lead to 
an asynchronous, non-oscillating steady state dominated 
by oscillators in state 1. For finite populations, fluctua¬ 
tions may also lead the population to stochastically jump 
from one stable state (oscillations) to the other (fixed 
point) (Figure [HI upper right). As k is further increased, 
the unstable fixed point at Pi Ri 1/3 becomes stable (su¬ 
percritical Hopf), leading to a small region of bistability 
between two non-oscillating states. Finally, at A: r; 3.4, 
the lower branch of the solution disappears and the pop¬ 
ulation settles into a fixed point on the upper branch of 
the solution curve. Similar behavior occurs for other y 
values in the approximate range 0.8 < y < 1 ; for smaller 
y, the Hopf bifurcation occurs prior to the emergence of 
the upper solution branch, leading to a larger region of 
bistability between asynchronous fixed points. 


1. Bistability in Other Classes of Oscillators 

When oscillators in any state can divide (i.e. Case 1), 
our linear stability analysis shows that similar behavior 
exists for a class of coupling functions, as long as the 
derivatives follow certain restrictions (see Section III Cl) . 
Unfortunately, because the rotational symmetry of the 
model is broken when only one state can divide (i.e. Case 
2 ), we are not able to provide similar analytical evidence 
that the bistable behavior discussed in Section Hi PI le.g. 
Figure El) occurs for other choices of coupling function. 
However, it is straightforward to perform numerical lin¬ 
ear stability analysis for any particular model. 

To confirm that these findings are not specific to the 
chosen form of the coupling function, we performed lin¬ 
ear stability analysis for F^ = at a = 3.65 > ac,o- 
As shown in Figure [71 we see similar dynamics for this 
choice of coupling. Specifically, for small k, we see stable 
oscillations. As k increases, we see a region of bistability 
between synchronous oscillations and an asynchronous 
fixed point, followed by a region of bistability between 
two asynchronous fixed points. Finally, at A: Ri 3, the 
lower branch of the solution disappears and the popula¬ 


tion settles into a fixed point on the other branch of the 
solution curve. Numerically, we find that similar bista¬ 
bilities also occur for other coupling functions, includ¬ 
ing multiple examples of the form F^ = 1 -|- aP" with 
1 < n < 4 (results not shown). As before, these coupling 
functions lead to a supercritical Hopf bifurcation in the 
absence of growth. 

Our analysis of discrete phase oscillators suggests that 
population growth can induce bistability between stable, 
asynchronous fixed points and stable, synchronous os¬ 
cillations when the growth is strongly phase-dependent. 
While the exact mathematical correspondence between 
these discrete phase models and classic continuous phase 
models, such as the Kuramoto model, has yet to be rigor¬ 
ously established, one might expect similar bistability to 
occur in continuous phase models as well. To explore this 
question, we performed numerical simulations of pop¬ 
ulations of Kuramoto oscillators [l|, whose (continous) 
phases evolve according to 

=u;i +- (j)i), (26) 

3 

uji is the intrinsic frequency of oscillator i, AT is a coupling 
parameter, and the sum runs over all oscillators in the 
population. We take uii = oj for all oscillators; for identi¬ 
cal oscillators, a synchronous state exists for all K > 0. 
To incorporate phase dependent population growth, at 
each time step we allow oscillators whose phase falls in 
a given range, i/o < < '/’i, to reproduce with prob¬ 

ability per unit time k. Following division, the daugh¬ 
ter oscillator has a phase that is chosen to fall in the 
range (j)o < (f ^ 4’i with uniform probability. To circum¬ 
vent numerical limitations due to exponentially growing 
populations, we start with Nq = 500 oscillators and al¬ 
low them to grow to a total size of A^ = 10®; when N 
reaches this maximum value, we reset the population 
size to A^ = 1000 while preserving the phase distribu¬ 
tion of oscillators. In practice, we preserve the approxi¬ 
mate phase distribution by binning oscillator phases into 
a total of M bins over the range 0 < ^ < 27r. Prior to 
resetting the population size, we calculate the fraction 
of oscillators in each bin and choose the phases in the 
smaller population so that the distribution is conserved. 
We find that similar behavior is observed as long as M 
is sufficiently large. In what follows, we choose M = 10, 
(po = 7r/3, cpi = 27r/3, and we set K = 1. 

As with the discrete oscillators, numerical simulations 
suggest that population growth decreases the synchrony 
of the population, eventually leading to a discontinuous 
transition to asynchronous behavior and a region of bista¬ 
bility between synchronous and asynchronous oscillations 
(Figure El). Because the simulations involve finite popu¬ 
lations, it is common to see trajectories that stochasti¬ 
cally switch between stable oscillations and asynchronous 
behavior (Figure El inset). While a full analysis of con¬ 
tinuous models is an exciting avenue for future work, we 
stress that our goal here is simply to provide evidence 
that growth-induced bistability can occur in continuous 
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phase models. Further analysis along these lines is nec¬ 
essary and would be welcome, but it is beyond the scope 
of the current work. 


III. EXCITABLE ELEMENTS 

A. Model for Growing Populations of Excitable 
Elements 

The results of Section [H] raise the question of whether 
population growth might have similar effects on popu¬ 
lations of coupled excitable elements. To explore this 
question in a simple context, we will model an excitable 
system as a discrete m-state system comprised of a qui¬ 
escence state ( 0 ), an excited state ( 1 ), and a finite set of 
refractory states (2,3, ...,m— 1) using the model intro¬ 
duced in [U. For each state i, the discretized phase is 
9 = Here we study a simple four state system, 

which contains the minimum number of refractory states 
(2) required for stable synchronization [s^. 

This discrete time model involves deterministic transi¬ 
tions between states 1 —>■ 2 and 2 —)> 3 and probabilistic 
transitions from the quiescent state 0 to the activated 
state 1. Coupling between elements is achieved by allow¬ 
ing the transition from state 0 to state 1 to depend on 
the states of the neighboring systems and a parameter a 
that measures the strength of the inter-element coupling. 
The last transition, from state 3 back to the quiescent 
state 0 , is also probabilistic and occurs with probability 
which is the same for all excitable elements [2^. A 
thorough analysis of this model [ 2 ^ reveals that, in the 
thermodynamic limit, the population transitions from an 
absorbing state to an active state at a critical value of 
a. More interestingly, for certain choices of p.y, it can 
undergo synchronous oscillations within a range of cou¬ 
pling strengths a\ < a < and there exists a region of 
bistability between synchronous oscillations and an asyn¬ 
chronous fixed point for CT 2 < ct < 173 , with cti < U 2 < 173 . 

In what follows, we study this model in the case where 
individual elements are allowed to reproduce, producing 
daughter elements in potentially different states. As with 
the oscillators, we focus on two cases: all oscillators are 
equally likely to divide (case 1 ), and division occurs only 
for oscillators in one particular state (case 2 ). As before, 
the daughter element will be in the same state as the 
mother with probability y, and will be randomly assigned 
to each of the other states with probability (1 —x)/3. The 
discrete nature of the model makes it amenable to rapid 
numerical simulations [ 2 ^ , which we explore in the next 
section. 


B. Numerical Simulations 

We performed discrete time simulations starting from 
Nq = 10^ globally coupled excitable elements for t = 500 
total time steps. To avoid limits due to computational 


memory, we reduce the total population size by a factor 
of 5 (while preserving the state distribution) when the 
number of units reaches 5 x 10"^. As in we take the 
probability of exciting a quiescent state to be 




Nt(l) 


(27) 


where N{t) is the total number of units at time step t, 
Nt{l) is the number of elements in state 1 at time step 
t, each excited element can activate a neighbor in the 
quiescent state with probability a/N{t), and reflects 
the probability that the quiescent state is excited by at 
least one of its Nt{l) active neighbors. We choose initial 
conditions so that Tb( 0 ) = 0.8 and Po(l) = 0 . 2 , where 
Pi{j) is the fraction of oscillators in state j at time step 
i. For each choice of k and %, we continuously increased 
coupling strength a from 2 to 50 for each successive sim¬ 
ulation, using the steady state from the previous value 
of a as initial conditions for the next value. We then 
repeated the simulations starting from cr = 50 and de¬ 
creasing a to 0 ; when simulations reached a steady state, 
we calculated the order parameter g, which is the same 
as that used for oscillators (see Section HB. 

In the case where all states can divide, the most salient 
effect of population growth is to decrease the range of cr 
over which oscillations are stable (Figure IH]). For ex¬ 
ample, when p^ = 0.94, 2^ showed that the model 
undergoes a synchronizing transition as a is increased 
from zero. Further increase of cr leads to discontinuous 
re-entrant transition that includes a region of bistabil¬ 
ity between synchronous and asynchronous states (black 
curves, main panel). When including population growth, 
we find that as k is increased from zero, the size of both 
the synchronized and bistable regions shrink (Figure IHl 
blue curves); as k is further increased, the bistable region 
eventually disappears (Figure [9l red curves) and, eventu¬ 
ally, the entire synchronization region is lost (not shown). 
We observe a similar decrease in the size of the oscillatory 
regime for x = 0 (top panels) and x = 1 (bottom pan¬ 
els) and for smaller values of p^ where bistability does 
not exist, even in the absence of growth. In all cases, 
increasing x slightly counteracts the effect of population 
growth; that is, larger x leads to a slightly larger region 
of synchrony and/or bistability. 

We find qualitatively similar behavior when division is 
restricted to only one state, such as state 1 iFigure ITOl). 
Interestingly, however, we find that the effect of increas¬ 
ing X will depend on which state is chosen for division. 
Specifically, when division is restricted to states 0, 2 or 3, 
increasing x will lead to an increase in the size of the ac¬ 
tive region iFigure fTUl upper insets), as in the case where 
all states can divide (Figure H]). By contrast, when di¬ 
vision is restricted to state 1 , increasing x will lead to a 
decrease in the size of the active region (Figure [TUI main 
panels and lower insets). 

For all simulations, we also calculated 1) the frequency 
of macroscopic oscillations in the active regions and 2 ) 
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FIG. 9. Order parameter q vs coupling strength a when all 
states divide for (a) x — 0 and (b) x = 1- Error bars are 
from standard deviations over 10 rnns. Black dashed, blue 
dash dot and red dotted line represent growth rate fc = 0, 
0.005 and 0.01, respectively, at pj = 0.94. For smaller values 
of Pj, the size of the active region shrinks and the bistable 
area disappears, even in the absence of growth (npper insets 
in (a) and (b): pj = 0.84). If pj is fnrther reduced, stable 
oscillations will eventually disappear (lower insets in (a) and 
(b): p^ = 0.74). 


FIG. 10. Order parameter q vs coupling strength a when 
division is restricted to one state. Main panels: Division is 
restricted to units in state 1: (a) x = 0 and (b) x = 1- Black 
dashed, blue dash dot and red dotted lines represent growth 
rate fc = 0, 0.01 and 0.07, respectively at p-y = 0.94. When 
increasing x from 0 to 1, the size of the active region shrinks 
(lower insets of both panels, k — 0.07). For comparison, upper 
insets show effect of increasing x when division is independent 
of state (upper insets in (a) and (b), k = 0.01). 


the maximum growth rate-which we refer to as the crit¬ 
ical growth rate-that still allows for synchronous oscil¬ 
lations. Figure ma) shows that population growth in¬ 
creases the frequency of oscillations, regardless of which 
state is chosen for division. Figure ITIT b) illustrates that 
as X increases, the critical growth rate increases when 
division is independent of state or is restricted to states 
0, 2, or 3; the opposite trend exists when division is re¬ 
stricted to state 1 (upper red triangles). 


C. Mean Field and Linear Stability Analysis 


To gain a systematic picture of these results, we follow 
the approach in to develop a mean field approxima¬ 
tion and derive full phase diagrams for these excitable 
systems. Specifically, by assuming that probability of 
exciting a quiescent state via one of its N{t) — 1 neigh¬ 
bors is aPt{l)/{N{t) — 1), the probability of excitation 
via at least one excited neighbor is 


PmW = 1 


\ m) -1 


(JV(t)-l) 


(28) 


In the thermodynamic limit, becomes 1 — In 

the case where all states divide, the system is governed 
by the difference equations 


Nt+iiO) = p^Nti3) + e-‘^^*(i)iVi(0) 

+ HxNtiO) + i^(fVt(l) + Nt{2) + Nt{3))) 

Nt+iil) = (l-e-"^‘(i))iVt(0) 

+ k{xNt{l) + + Nti2) + Ntim 

Nt+i{2) = Nt{l) 

+ k{xNt{2) + i^(fVt(0) + Nt{l) + Nt{3))) 
Nt+i{3) = Nt{2) + {1 - p^)Nt{3) 

+ HxPtii) + + fVt(l) + fVt(2))), 

(29) 


and the total number of oscillators follows N{t+1) = (H- 
k)N{t). If all equations are divided by the total number 
of oscillators, we arrive at equations for the probability 
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FIG. 11. Frequency and critical growth rate vs X- Upper 
figure: circle, square and up triangle represent x = 0, 0.5 and 
1, respectively. Blue color shows the case in which all states 
divide and red color means only state 1 divides, = 0.9 and 
CT = 8. Lower figure: black circle, blue square, red up triangle, 
magenta down triangle and cyan diamond represent hve cases 
(all states divide, only state 0 divides...), respectively. Inset: 
zoom in black circle (all states divide case) 


Similarly, if division is restricted to only one state, such 
as state 1, we have 




Pj 


m) 


1 + kPtil) 

fc(l - x) 

■ 3(l + fcPt(l)) 


1 + kPt{l) 

m) 


PtiO) 


(34) 


1 — p 


l + A:Pt(l)^‘^^^^ 3(l + fcPt(l))^‘^^^ 


(35) 

(36) 


^^*+i(3) = 


Ptm 


1 + kPtil) 

fc(i - x) 

■ 3(l + fcPi(l)) 


1 -P7 

1 + kPtil) 
Pti^)- 


Pt{^) 


(37) 


It is easy to derive similar equations when division is 
restricted to one of the other states. In all cases, non¬ 
trivial steady states P*{i) occur when Pt+i{i) = Pt{i) 
{i = 0,1, 2, 3). Since the probabilities are normalized to 
1, we can omit equations (|3ni) and (|Mll and reduce each 
set of equations by one. 

To study the stability of nontrivial solutions, we lin¬ 
earize the equations near each solution; the correspond¬ 
ing Jacobian matrix J when all states can divide is given 

by 


\ + k 



(38) 


Pt{i) to be in state i at time t, 

P..r(0) = ^P.(3) + ^P.(0) 

^ (xPt(0) + i^(l-Pt(0))) 


(30) 


1 -I- A: 


where Jn is the Jacobian for non-dividing popula¬ 
tions HSl 


an - 1 - 1 

Jn “ ( 1 b 0 

0 1 1 — 7*7 


(39) 


p* 


1 _ p-<rPt(l) 

P.pr(l) = ^^P.(0) 

+ Y^ixPta) + ^{^-Ptm 

P.pr(2) = ^P*(l) 

+ Y^ixPt{2) + ^^{i-Ptm)) 

P*Pr(3) = ^P.(2) + ^P.(3) 

+ Y^ixPt{3) + ^^a-Ptim- 


(31) 


(32) 


(33) 


I is identity matrix, and an = ae~’^P {1 — Pi — P 2 — P 3 ) + 
— 1. At steady state, P* = P 2 = P 7 P 3 . We note 
that even when x = 1 the excitable elements repro¬ 
duce to form identical daughter cells, the model does not 
reduce to the corresponding non-growing model. Hence, 
unlike the oscillator model, growth will modify the dy¬ 
namics of these excitable systems even in the case when 
all states can divide and x = 1 ■ For cases where division 
is restricted to one state, the Jacobian can be readily cal¬ 
culated numerically, but it cannot be simply written in 
terms of Jat. 

As in the case of continuous time systems, the eigen¬ 
values provide information about the stability of each 
fixed point. Specifically, bifurcations between stable and 
unstable fixed points occur when |A| = 1, allowing us 
to dilineate phase boundaries separating oscillatory and 
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FIG. 12. Phase diagrams when all oscillator states can divide. 
Upper panel: Phase diagram vs a when all states divide 
and k = 0.005, x = 1- Curves indicate nature of bifurcation 
(blue, supercritical Hopf; red, subcritical Hopf; black, global 
saddle node of limit cycles). Regions of bistability between 
synchronous and asynchronous states are indicated. Crosses 
and circles indicate results from numerical simulations (A^o = 
10"^). Black line shows numerical ctcs from mean field analysis. 
Cross and circle represent simulation results. Dash dot means 
non-growth for comparison. Lower panels: phase boundaries 
(excluding bistable regions) when all states divide, x = Oj 
fc = 0 (black circle), k = 0.01 (blue square), k = 0.02 (red 
up triangle), k — 0.03 (magenta down triangle). Lower right 
inset, k = 0.02 for different values of X- Black circle, blue 
square, red up triangle, magenta down triangle represent x = 
0, |, I and 1, respectively. 


non-oscillatory regimes. Specifically, we are interested in 
the location of Neimark-Sacker (NS) bifurcations, which 
are the discrete time analog of a Hopf bifurcation. As 
with Hopf bifurcations, the nature of the transition is 
given by the sign of the first Lyapunov coefficient: the 
NS bifurcation is supercritical if Zi < 0 and subcritical if 
li > 0. Following standard bifurcation theory (see, for 
example, M) , we calculate li as 


^i = ^^e{A (^s,C(f,r,r)) 

+ 2(s,.B(r,(/-J)-iH(f,^); 

+ (s,H(F, (A^J - J)-iH(r,r)))] } , 


(40) 


where brackets represent the standard complex inner 



FIG. 13. Phase diagrams when division is restricted to one 
state. Main panel: phase diagram for excitable elements with 
growth (solid lines; x = 0, A: = 0.04, only state 1 can divide) 
and without growth (dashed lines). Curves indicate nature 
of bifurcation (blue, supercritical Hopf; red, subcritical Hopf; 
black, global saddle node of limit cycles). Regions of bista¬ 
bility between synchronous and asynchronous states are in¬ 
dicated. Crosses and circles indicate results from numerical 
simulations {No = 10^). Right panel, phase boundaries (ex¬ 
cluding bistable regions) for x = 0 and X = 1 when only state 
1 (top) or only state 2 (bottom) divides; k = 0.07 in both 
cases. 


product, r is the eigenvector of the Jacobian matrix with 
corresponding eigenvalue A, s is the eigenvector of 
with eigenvalue A*, and normalization is chosen so that 
(r, r) = 1 and (s, r) = = 1- Furthermore, 

B{x, y) and C{x, y, z) are vector-valued multi-linear func¬ 
tions 


3 

B^{x,y) = 


j,k=l 


d^F,{Pt) 

dPt{j)dPt{k) 


xjyk, 

P* ,rTc 


C^{x,y,^= Y 


d^P{Pt) 


dPt{j)dPt{k)dPt{l) 


xjykzi, 


P* ,<^c 


where x = (cci, 0 : 2 , 2 : 3 )^, y = (yi,j/ 2 ^)^ and z = 
(zi,Z 2 , zs)'^ are arbitrary vectors. As in [ 2 ^, when li > 0 
(indicating a subcritical bifurcation and the correspond¬ 
ing bistability), we supplement the above calculations 
with simulations of the mean field equations to deter¬ 
mine the location of the bistable regimes, which are not 
fully determined by linear stability properties. 

Figure [ 12 ] shows the resulting phase diagrams for the 
case where all states divide, and Figure [13] shows the 
phase diagram when division is restricted to state 1. Pop¬ 
ulation growth has several obvious effects on the excitable 
systems. First, growth can shift the location of the fixed 
points, but unlike the oscillator case, we do not find ev¬ 
idence of new fixed points. Consistent with our simula¬ 
tions, growth shifts the phase boundaries to higher values 

























14 


of leading to a smaller region of synchronized oscil¬ 
lations and a significantly reduced region of bistability; 
as k is further increased, the oscillatory region is even¬ 
tually eliminated (see Figure [T^l bottom panel). The 
phase diagrams indicate that adding growth can have 
significant effects on the dynamics, depending on the val¬ 
ues of a and p-y. Consider, for example, the main panel 
of Figure [T31 Systems originally undergoing oscillations 
can enter a bistable state (e.g. {(7,p^) = 18,0.98)) or a 
non-oscillating active state (e.g. {(J,pj) = 8,0.82)) when 
growth is increased to fc = 0.04. Qualitatively similar 
behavior is observed when division is restricted to one of 
the other states, or when all states can divide fFigurelT^. 

Interestingly, we also find that the effect of x on the 
phase diagram will depend on which state is chosen for di¬ 
vision (Figure[T21 bottom inset; Figure[T21 right panels). 
For example, increasing x at a fixed value of k will raise 
the phase boundary to higher p^ when state 1 divides, 
but will lower the boundary when states 0, 2, or 3 divide 
(Figure [T31 right panels) or when all states can divide 
(Figure [121 bottom inset). When only the excited state 
(1) divides, self-similarity between mother and daughter 
oscillators decreases the area of phase space over which 
synchrony can occur, while such self-similarity increases 
the size of the synchronized region when any of the non- 
excited states divide. 


IV. DISCUSSION 

We have shown that population growth, and the corre¬ 
sponding redistribution of oscillator phases, can induce 
a wide range of new dynamic behaviors in systems of 
coupled oscillators and excitable elements. In particu¬ 
lar, when growth is independent of oscillator phase, in¬ 
creasing growth rate leads to an increase in oscillation 
frequency and a decrease in phase synchrony, eventually 
culminating in a transition to a non-oscillating steady 
state. Interestingly, the growth-induced transition can be 
subcritical, even when the non-growing model exhibits a 
supercritical bifurcation; in that case, one sees a bistable 
region with coexisting synchronous and asynchronous 
states, depending on the derivatives of the coupling func¬ 
tion. When division is strongly state dependent, growth 
can again lead to extinction of oscillations, but in certain 
parameter regimes, one finds new asynchronous states 
with unequal phase distributions, bistability between two 
asynchronous states or between asynchronous and oscil¬ 
latory states, or modulation of the bulk oscillation fre¬ 
quency. In excitable systems, which may include bistable 
behavior even in the absence of growth, the most salient 
effects of growth are to shift the phase boundaries sepa¬ 
rating active, oscillatory, and bistable regimes while also 
increasing the frequency of super-threshold oscillations. 
In practice, the shifting phase boundaries can lead to a 


range of different dynamical effects, many of which mirror 
the behaviors seen in oscillators. For example, systems 
originally undergoing oscillations can enter a bistable 
state or a non-oscillating active state when growth is in¬ 
creased. 

Several recent studies have focused on related ideas. 
First, a number of studies have examined cell-density de¬ 
pendent synchronization (see, for example, [Il,[il). In 
these systems, synchronization of intracellular dynam¬ 
ics can be modified at high cell densities as a result 
of biochemical communication known as quorum sens¬ 
ing. These studies do not directly explore the effects 
of population growth rate, but instead treat cell density 
as the relevant control parameter that governs not only 
intracellular coupling, but also the dynamics of individ¬ 
ual cells. In addition, the authors of develop a de¬ 
tailed biochemical model of circadian clocks in growing 
cyanobacteria populations. They show that otherwise 
stable oscillations-specifically, those driven by phospho¬ 
rylation cycles-are destabilized at high growth rates via a 
supercritical Hopf bifurcation for the chosen range of pa¬ 
rameters. Therefore, in fast growing populations, addi¬ 
tional stabilizing mechanisms (transcription-translation 
cycles) are required to preserve integrity of the oscilla¬ 
tions. It would be interesting to see if subcritical bifur¬ 
cations, including the bistability observed in our models, 
occur in different parameter regimes. In bacteria popu¬ 
lations, an elegant series of studies has recently laid the 
groundwork for synthetic sensors and logical program¬ 
ming in living systems based on tunable spatiotemporal 
oscillations in growing populations [l^ . [l7j |. However, 
the effects of growth rate on these dynamics are not ad¬ 
dressed in detail. Finally, recent theoretical work [s^ has 
demonstrated a rich collection of dynamical behaviors in 
small chains of coupled oscillators gradually increasing 
in number. By contrast, here we focus on large systems 
of oscillators and the corresponding phase transitions to 
macroscopic oscillations. To our knowledge, this work 
is the first to systematically address how coupling be¬ 
tween population growth and oscillator phase can effect 
synchronization. 

We have shown that population growth can dramati¬ 
cally influence synchronization phenomena and, in some 
cases, lead to entirely new dynamical states in popula¬ 
tions of coupled oscillators and excitable elements. While 
we focused on discrete phase models because of their rel¬ 
ative simplicity, we hope these results will motivate fu¬ 
ture explorations on the interplay between synchroniza¬ 
tion and population growth in additional models of os¬ 
cillators, excitable elements, and perhaps more general 
dynamical systems. Given the theoretical importance of 
self-synchronization in statistical physics and its ubiquity 
in biological systems, we believe that the potential effects 
of population growth on collective oscillations will prove 
to be an important and rich topic for future exploration. 
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